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ABSTRACT 

We introduce a novel Earth-like planet surface temperature model (ESTM) for habitabil¬ 
ity studies based on the spatial-temporal distribution of planetary surface temperatures. The 
ESTM adopts a surface Energy Balance Model complemented by: radiative-convective atmo¬ 
spheric column calculations, a set of physically-based parameterizations of meridional transport, 
and descriptions of surface and cloud properties more refined than in standard EBMs. The 
parameterization is valid for rotating terrestrial planets with shallow atmospheres and moder¬ 
ate values of axis obliquity (e < 45°). Comparison with a 3D model of atmospheric dynamics 
from the literature shows that the equator-to-pole temperature differences predicted by the two 
models agree within « 5 K when the rotation rate, insolation, surface pressure and planet ra¬ 
dius are varied in the intervals 0.5 < D/De ^ 2, 0.75 < S/So ^ 1.25, 0.3 < p/(lbar) < 10, 
and 0.5 < R/R® % 2, respectively. The ESTM has an extremely low computational cost and 
can be used when the planetary parameters are scarcely known (as for most exoplanets) and/or 
whenever many runs for different parameter configurations are needed. Model simulations of 
a test-case exoplanet (Kepler-62e) indicate that an uncertainty in surface pressure within the 
range expected for terrestrial planets may impact the mean temperature by ^ 60 K. Within the 
limits of validity of the ESTM, the impact of surface pressure is larger than that predicted by 
uncertainties in rotation rate, axis obliquity, and ocean fractions. We discuss the possibility of 
performing a statistical ranking of planetary habitability taking advantage of the flexibility of 
the ESTM. 


Subject headings: planetary systems - astrobiology 

1. Introduction 


The large amount of exoplanet data collected 
with the doppler and transit methods (e.g. [Mayor 


et al.|[20l4 Batalha||2013 and refs, therein) indi¬ 


cate that Earth-size planets are intrinsically more 
frequent than giant ones, in spite of the fact that 
they are more difficult to detect. Small planets are 
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found in a relatively broad range of metallicities 
(Buchhave et al. and, at variance with giant 

planets, their detection rate drops slowly with de¬ 
creasing metallicity (Wang & Fischer 20151. These 
observational results indicate that Earth-like plan¬ 


ets are quite common around other stars (e.g. Farr 


et al. 20141 and are expected to be detected in 


large numbers in the future. Their potential sim¬ 
ilarity to the Earth makes them primary targets 
in the quest for habitable environments ouside the 
Solar System. Unfortunately, small planets are 
quite difficult to characterize with experimental 
methods and a significant effort of modelization is 
required to cast light on their properties. The aim 
of the present work is to model the surface tem¬ 
perature of these planets as a contribution to the 
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study of their surface habitability. The capability 
of an environment to host life depends on many 
factors, such as the presence of liquid water, nu¬ 
trients, energy sources, and shielding from cosmic 
ionizing radiation (e.g. Seager|[20l3 Guedel et al. 
20141. A knowledge of the surface temperature 


is essential to apply the liquid water criterion of 
habitability and can also be used to assess the po¬ 
tential presence of different life forms according to 
other types of temperature-dependent biological 
criteria (e.g. Clarke|[2014[ ). Here we are interested 
in modeling the latitudinal and seasonal variations 
of surface temperature, T[ip, t), as a tool to calcu¬ 
late temperature-dependent indices of fractional 
habitability (e.g. Spiegel et al.||2008 ). 

Modeling T{ip,t) is a difficult task since many 
of the physical and chemical quantities that gov¬ 
ern the exoplanet surface properties are currently 
not measurable. A way to cope with this prob¬ 
lem is to treat the unknown quantities as free pa¬ 
rameters and use fast climate calculations to ex¬ 
plore how variations of such parameters affect the 
surface temperature. General Circulation Models 
(GCMs) are not suited for this type of exploratory 
work since they require large amounts of compu¬ 
tational resources for each single run as well as a 
detailed knowledge of many planetary character¬ 
istics. Two types of fast climate tools are com¬ 
monly employed in studies of planetary habitabil¬ 
ity: single atmospheric column calculations and 
energy balance models. Atmospheric column cal¬ 
culations treat in detail the physics of vertical en¬ 
ergy transport, taking into account the influence 
of atmospheric composition on the radiative trans¬ 
fer (e.g., Kasting|1988 |. This is the type of climate 
tool that is commonly employed in studies of the 


“habitable zone” (e.g., Kasting 1988 Kasting et al. 


1993[|von Paris et al.|2013 Kopparapu et al.|2013 


20141. Energy balance models (EBMs) calculate 


the zonal and seasonal energy budget of a planet 
using a heat diffusion formalism to describe the 
horizontal transport and simple analytical func¬ 
tions of the surface temperature to describe the 
vertical transport (e.g.. North et al.||1981 ). EBMs 
have been employed to address the climate im¬ 
pact induced by variations of several planet pa¬ 
rameters, such as axis obliquity, rotation period 


and stellar insolation (Spiegel et al. 2008 2009 


Dressing et al.|[T010[ [Spiegel et al. 2010 Eorgan 
2012, 20141. By feeding classic EBMs with multi¬ 


parameter functions extracted from atmospheric 
column calculations one can obtain an upgraded 
type of EBM that takes into account the physics of 
vertical transport (Williams & Kasting 19971. Fol¬ 
lowing a similar approach, in a previous paper we 
investigated the impact of surface pressure on the 
habitability of Earth-like planets by incorporat¬ 
ing a physical treatment of the thermal radiation 
and a simple scaling law for the meridional trans¬ 
port ( Vladilo et al.|[20T3 hereafter Paper 1). Here 
we include the transport of the short-wavelength 
radiation and we present a physically-based treat¬ 
ment of the meridional transport tested with 3D 
experiments. In this way we build up an Earth¬ 
like planet Surface Temperature Model (ESTM) 
in which a variety of unknown planetary proper¬ 
ties can be treated as free parameters for a fast 
calculation of the surface habitability. The ESTM 
is presented in the next section. In Section 3 we 
describe the model calibration and validation. Ex¬ 
amples of model applications are presented in Sec¬ 
tion 4 and the conclusions summarized in Section 
5. 


2. The model 

The ESTM consists of a set of climate tools 
and algorithms that interchange physical quanti¬ 
ties expressed in parametric form. The core of the 
model is a zonal and seasonal EBM fed by multi¬ 
parameter physical quantities. The parameteri¬ 
zation is obtained using physically-based climate 
tools that deal with the meridional and vertical 
energy transport. The relationship between these 
ingredients is shown in the scheme of Fig. In 
the following we present the components of the 
model, starting from the description of the EBM. 

In zonal EBMs the surface of the planet is di¬ 
vided into zones delimited by latitude circles. The 
surface quantities of interest are averaged in each 
zone over one rotation period. In this way, the 
spatial dependence is determined by a single coor¬ 
dinate, the latitude ^p. Since the temporal depen¬ 
dence is “smoothed” over one rotation period, the 
time, t, represents the seasonal evolution during 
the orbital period. The thermal state is described 
by a single temperature, T = representa¬ 

tive of the surface conditions. By assuming that 
the heating and cooling rates are balanced in each 
zone, one obtains an energy balance equation that 
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Fig. 1.— Scheme of the Earth-like planet Surface Temperature Model (ESTM). The zonal and seasonal 
EBM (central box) is fed by physical quantities (circled symbols) described in §2. At variance with classic 
EBMs, the physical quantities are multi-parameter functions modelled with the aid of atmospheric column 
calculations (I and A) and 3D climate experiments (D). 


is used to calculate The most common 


form of EBM equation ( 

North et al. 

1981 

Williams 

& Kastingl 19971 ISpiege 

letal.|2008 

Pierrehumbert 


2010 Gilmore||2014) is 


dT d 

dt dx 


+ / = 5(1-A), (1) 


where x = sin tp and all terms are normalized per 
unit area. The first term of this equation repre¬ 
sents the zonal heat storage and describes the tem¬ 
poral evolution of each zone; C is the zonal heat 
capacity per unit area (North et al. 1981). The 
second term represents the amount of heat per 
unit time and unit area leaving each zone along 
the meridional direction (North et al. 1981 Eq. 
(21)). It is called the “diffusion term” because the 
coefficient D is defined on the basis of the analogy 
with heat diffusion, i.e. 


$ = 

0(p 


( 2 ) 


where 27ri?^d> cos (p is the net rate of energy trans¬ 
port^ across a circle of constant latitude and R is 


^ With the adopted definitions of D and <I> it is easy 


the planet radius (see Pierrehumbert 2010). The 
term I represents the thermal radiation emitted 
by the zone, also called Outgoing Longwave Radi¬ 
ation (OLR). The right side of the equation repre¬ 
sents the fraction of stellar photons that heat the 
surface of the zone; S is the incoming stellar radi¬ 
ation and A the planetary albedo at the top of the 
atmosphere. All coefficients of the equation de¬ 
pend, in general, on both time and latitude, either 
directly or indirectly, through their dependence on 
T. 


In classic EBMs the coefficients D, I and A are 
expressed in a very simplified form. As an ex¬ 
ample, D is often treated as a constant, in spite 
of the fact that the meridional transport is influ¬ 
enced by planetary quantities that do not appear 
in the formulation ([^. The OLR and albedo are 
modelled as simple analytical functions, / = I{T) 
and A = A(T), while they should depend not 
only on T, but also on other physical/chemical 
quantities that influence the vertical transport. 


to show that ^(<I>cos(/j) = — “ 

— ^ 1^(1 — represents the latitudinal transport 

per unit area (see|North et al.||1981[ Eq. 21). 
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This simplified formulation oi D, I and A pre¬ 
vents important planetary properties to appear 
in the energy balance equation 0. To obtain 
a physically-based parameterization we describe 
the vertical transport using single-column atmo¬ 
spheric calculations and the meridional transport 
using algorithms tested with 3D climate experi¬ 
ments. Thanks to this type of parameterizatiorQ 
the ESTM features a dependence on surface pres¬ 
sure, p, gravitational acceleration, g, planet ra¬ 
dius, R, rotation rate, D, surface albedo, as, stellar 
zenith distance, Z, atmospheric chemical composi¬ 
tion, and mean radiative properties of the clouds. 
By running the simulations described in Appendix 
A, the ESTM generates a “snapshot” of the sur¬ 
face temperature T{(p, t) in a very short computing 
time, for any combination of planetary parameters 
that yield a stationary solution of Eq. Q. We now 
describe the parameterization of the model. 


2.1. The meridional transport 

The heat diffusion analogy 0 guarantees the 
existence of physical solutions and contributes to 
the high computational efficiency of EBMs. In or¬ 
der to keep these advantages and at the same time 
introducing a more realistic treatment of the lati¬ 
tudinal transport, here we derive and D in terms 
of planet properties relevant to the physics of the 
horizontal transport. To keep the problem simple 
we focus on the atmospheric transport (the ocean 
transport is discussed below in { 2.1.3| ) . The atmo¬ 
spheric flux can be derived applying basic equa¬ 
tions of fluid dynamics to the energy content of a 
parcel of atmospheric gas. The energy budget of 
the parcel is expressed in terms of the moist static 
energy (MSE) per unit mass, 


m = CpT -I- LyVv -f gz 


( 3 ) 


where the terms CpT, LyVy and gz measure the 
sensible heat, the latent heat and the potential en¬ 
ergy content of the parcel at height z, respectively; 
Ly is the latent heat of the phase transition be¬ 
tween the vapor and the condensed phase; Vy the 
mass mixing ratio of the vapor over dry compo¬ 
nent; g is the surface gravity acceleration. The 
MSE and the velocity of the parcel are a function 


of time, t, longitude. A, latitude, p, and height, 
z. The latitudinal transport is obtained by inte¬ 
grating the fluid equations in longitude and ver¬ 
tically, the height z being replaced by the pres¬ 
sure coordinate, p = p{z). Starting from a sim¬ 
plified mass continuity relation valid for the case 
in which condensation takes away a minimal at¬ 
mospheric mass ( Pierrehumbert|2010 §9.2.1), one 
obtains the mean zonal flux 


= ^ T" dX r vm—= ^-vm (4) 
RJo Jo 9 Rg 

where p is the surface atmospheric pressure and v 
the meridional velocity component of the parcel. 
The second equality of this expression is valid for 
a shallow atmosphere, where g can be considered 
constant, as in the case of the Earth. 

To proceed further, we assume that 0 is valid 
when the physical quantities are averaged over one 
rotation period, since this is the approach used in 
EBMs. In this case, the time t represents the sea¬ 
sonal (rather than instantaneous) evolution of the 
system: variability on time scales shorter than one 
planetary day are averaged out. At this point we 
split the problem in two parts. First we derive a 
relation for $ and D valid for the extratropical 
transport regime. Then we introduce a formalism 
to empirically improve the treatment of the trans¬ 
port inside the Hadley cells. 


2.1.1. Transport in the extratropical region 

We consider an ideal planet with constant inso¬ 
lation and null axis obliquity, in such a way that 
we can neglect the (seasonal) dependence on t. We 
restrict our problem to the atmospheric circulation 
typical of fast-rotating terrestrial-type planets, i.e. 
with latitudinal transport dominated by eddies in 
the baroclinic zone. A commonly adopted formal¬ 
ism used to treat the eddies consists in dividing 
the variables of interest into a mean component 
and a perturbation from the mear0 , representa¬ 
tive of the eddies. By indicating the mean with an 
overbar and perturbations with a prime, we have 
for instance v = v + v' and m = m + m'. It is 
easy to show that vm = vfn v'm'. When the 
eddies transport dominates, the term vm can be 


A simpler parameterization, not tested with 3D climate 
calculations, was adopted by |Williams &: Kasting] | |1997^ 
and in Paper I. 


In general, the “mean” and the “perturbations” are re¬ 
ferred to time and or space variations. 
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neglected so that 


1 p- 


$ ~ - v'm' 

Rg 


and we obtain 


( 5 ) 


£» = -$ 



J_P 

R^g 



v'm' 


( 6 ) 


where dy = R dp is the infinitesimal meridional 
displacement. To calculate v'm' we consider the 
surface value of moist static energ 30 m = CpT + 
LyVy, from which we obtain 

v'm' = CpV'T' + Lyv'r'y . (7) 


We express the mean values of the perturbation 
products as 

( 8 ) 


v'T' = ks \v\ \T 


and 


v'r'^ = ki,\v'\\ry\, (9) 

where 11 means a root-mean square magnitud^ 
and ks and fcn are correlation coefficients (e.g. 


Barry et al. 2002 1 . At this point we need to 


quantify the perturbations of T and r„, i.e. of 
the quantities being mixed. In eddy diffusivity 
theories these perturbations can be written as a 
mixing length, fmix, times the spatial gradient of 
the quantity. We consider the gradient along the 
meridional coordinate y and we write 


and 


dT 

IT'I - — 

( 10 ) 

\' V\ - -^miX rv 

dy 

( 11 ) 


where T and are mean zonal quantities; since 
the mixing is driven by turbulence we assume that 
the mixing length is the same for sensitive and 
latent heat. To estimate dr^/dy we recall that 


dv Pv Pv qPv 

Tv = -=- 

Pdiy Pdry Pdry Pdry 


( 12 ) 


^The MSE is conserved under conditions of dry adiabatic as¬ 
cent and is approximately conserved in saturated adiabatic 
ascent. Therefore the MSE is, to some extent, independent 
of Results obtained by |Lapeyre Sz Held] < (200^ suggest 
that lower layer values of moist static energy are most ap¬ 
propriate for diffusive models of energy fluxes. 

^Root-mean square values must be introduced since the time 
mean of the linear perturbations is zero. 


where and pv are the molecular weight and 
pressure of the vapor, /Xdry and pdry the corre¬ 
sponding quantities of the dry air, q is the relative 
humidity and p* = p^(T') is the saturation vapor 
pressure. We assume constant relative humidity 
and we can write 


dry _ kdry dT\ _ q dpi dT 

dy \dT dy J pdry Pdry dT dy 


Combining the expressions from Q to (131 we ob- 
tain 


v'm' = —t 


■j'l ^ I 


khLi 


Pv 


q dpi 


Pdry Pdry 


dT 


and inserting this in ([^ we derive 


D ~ \v 

R^ g 


ksVp “t” kj^Ly 




q dpi 


Pdry Pdry dT 

At this point, we need an analytical expression 
for Among a large number of analyti¬ 

cal treatments of the baroclinic circulation (e.g., 
Green |1970[ |Stone|[T97^ [Gierasch fc Toon]|1973[ 


(14) 


.(15) 


Held 1999), here we adopt a formalism proposed 


by Barry et al. (2002) which gives the best agree¬ 
ment with GCM experiments. 


According to Barry et al. (2002 1 , the baroclinic 
zone works as a diabatic heat engine that obtains 
and dissipates energy in the process of transport¬ 
ing heat from a warm to a cold region. If we call Ty, 
and Ty the temperatures of the warm and cold re¬ 
gions, the maximum possible thermodynamic effi¬ 
ciency of the engine is 5T /T^,, where 5T = Tyy—Ty. 
The energy received by the atmosphere per unit 
time and unit mass, Q, represents the diabatic 
forcing of the engine. The rate of generation (and 
dissipation) of eddy kinetic energy per unit mass 
is given by 

'ST' 


e = r] 


7f^] Q 

-L 1/1 


(16) 


where t] is an efficiency factor representing the 
fraction of the generated kinetic energy used by 
heat transporting eddies. Assuming that the av¬ 
erage properties of the flow depend only on the 
length scale and the dissipation rate per unit 
mas^ dimensional arguments yield the velocity 
scaling law 


OC {ein 


■, 1/3 


(17) 


^If the eddies exist in an inertial range, the average prop¬ 
erties of the flow will depend only on the dissipation rate 
and the length scale ( [Barry et al.|200^ . 
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As far as the mixing length is concerned, the 
Rhines scale is adopted 

1/2 

(18) 


2\v'\ 


P 

where /3 = df /dy is the gradient of the Coriolis 
parameter, / = 2r2sin(/?, and the angular rota¬ 
tion rate of the planet. The study of [Barry et al.| 
(20021 suggests that, among other types of length 
scales considered in literature, the Rhines scale 
yields the best correlations in 3D atmospheric ex¬ 
periments. The adoption of the Rhines scale is 
also supported by a study of moist transport per¬ 


formed with GCM experiments (Frierson et al. 


20071. The Rhines scale must be calculated at the 
latitude of maximum kinetic energy, i.e. for 
P = (2Dcos(/3m)/7?. From the above expressions 
we obtain 

3/5 / ^ \ 




fl COS p-a 


Inserting this in (151 we obtain 


where 


D — Ddry (1 + A) 


Ddry = fcsCp7?^/^(cOS(pin) X 


9 


ST 

T,, 


Q 


3/5 


(19) 


( 20 ) 


( 21 ) 


is the dry component of the atmospheric eddies 
transport and 

khLv Pv q dpi 


A = 


( 22 ) 


ksCp //dry Pdry dT 

is the ratio of the moist over dry components. 

For the practical implementation of the analyt¬ 


ical expressions (20), (21) and (22) in the EBM 


code, we proceed as follows. The maximum ther¬ 
modynamic efficiency STjT^ is calculated by tak¬ 
ing T^ = T{pi) and Tc = T{p 2 ) where pi and 
P 2 are the borders of the mid-latitude region and 
overbars indicate zonal annual means. Follow¬ 


ing Barry et al. (2002), we adopt pi = 28° and 
P 2 = 68°, after testing that the model predic¬ 
tions are virtually unaffected by the exact choice 
of these value^U We estimate the diabatic forc¬ 
ing term (W/kg) as Q ~ {ASR}/(p/g), where 


r Also the GCM experiments by Barry et al. 


(2002 \ indicate 


that the results are not sensitive to the choice of ipi and 
V^2- 


{ASR} = {S'(l — A)} is the absorbed stellar ra¬ 
diation (W/m^) averaged over one orbital period 
in the latitude range (pi, P 2 ) and p/g the atmo¬ 
spheric columnar mass (kg/m^). We neglect the 
contribution of surface fluxes of sensible heat since 
they cannot be estimated in the framework of the 
EBM model. This approximation is not critical 
because these fluxes yield a negligible contribution 
to Q according to |Barry et al. (2002). Treating fcp, 
/cs, rj, and as constant^ we obtain from Eq. 


(21) a scaling law for the dry term of the transport 


§dry OC Cp i?-®/® ( H 0-4/® 




X 


3/5 
X 


(23) 


We estimate the temperature gradient of satu¬ 
rated vapor pressure as dpl/dT ~ Spl/ST, with 
^pI = \pl(Tw) -pliTc)]. Since fep, ks and are 
constants, we obtain from Eq. ( |2^ a scaling law 
for the ratio of the moist over dry components 


Snid 


Spl 


Cp /^dry Pdry 


ST 


(24) 


Finally, by applying Eq. (20) and the scaling laws 


(23) and (24) to a generic terrestrial planet and to 


the Earth, indicated by the subscript o, we obtain 


D 


§dry 

Sdry,o 


f T Aq * (Sjnd/^md,o) 
1 -l- Ao 


(25) 


With the above expressions we calculate D treat¬ 
ing R, O, p, g, Cp, //dry, Pdry, 9 as parameters 
that can vary from planet to planet, in spite of 
being constant in each planet. The ratio of moist 
over dry eddie transport of the Earth is set to 
Ao = 0.7 (e.g. Kaspi & Showman 2014). For 


the sake of self-consistency, we adopt the param¬ 
eters {St)o, {Tw)o, {ASR}o and (i5p*)o obtained 
from the Earth’s reference model. Since these 
parameters vary in the course of the simulation, 
we perform the calibration of the Earth model in 
two steps. First we calibrate the model excluding 
the ratiot]^ 5 t/(<5t)o, T^/T^^o, {ASR}/{ASR}„ 


^Numerical experiments performed with simplified GCMs 
suggest that the correlation coefficients /cl, kg and the ef¬ 
ficiency factor ?7 can be treated as constants with good ap¬ 
proximation ( [Barry et al. [200^ [Frierson et al.|2007^ . 

^Excluding these ratios is equivalent to setting them equal to 
unity in the Earth’s model, as they should be by definition. 
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and dpl/{6pl)o from the scaling laws of Eq. (251. 
Then we reintroduce these ratios in the scaling law 
adopting for {St)o, {Tw)o and {ASR}q the values 
(St), {Tyj), {ASR} and {5p%) obtained in the first 
step. The second step is repeated a few times, un¬ 
til convergence of the parameters {St)o, (Tiu)o and 
{ASR}(, and (dp*)o is achieved. 


2.1.2. Transport in the Hadley Cell 

The derivation performed above ignores the ex¬ 
istence of the Hadley Cells, since they do not con¬ 
tribute to the extratropical meridional transport. 
However, the Hadley circulation is extremely ef¬ 
ficient in smoothing temperature gradients inside 
the tropical region. This aspect cannot be com- 
pleteley ignored in our treatment, since our goal 
is to estimate the planet surface temperature dis¬ 
tribution. Unfortunately, the diffusion formalism 
of Eq. ([^ is inappropriate inside the Hadley Cells 
and the only way we have to improve the descrip¬ 
tion of the tropical temperature distribution is to 
correct the formalism with some empirical expres¬ 
sion. We summarize the approach that we follow 
to cope with this problem. 

The global pattern of atmospheric circulation 
is influenced, among other factors, by the sea¬ 
sonal variation of the zenith distance of the star. 
In the case of the Earth, a well known example 
of this type of influence is the seasonal shift of 
the Intertropical Convergence Zone (ITCZ), that 
moves to higher latitudes in the summer hemi¬ 
sphere. The ITCZ is, in practice, a tracer of the 
thermal equator at the center of the system of the 
two Hadley Cells, where we want to improve the 
uniformity of the temperature distribution. A way 
to do this is to enhance the transport coefficient D 
in correspondence with such thermal equator. To 
incorporate this feature in our model, we scale D 
according to mean diurnal value of t) = cos Z, 
where Z is the stellar zenith distance. In practice, 
we multiply Z? by a dimensionless modulating fac¬ 
tor, C,[ip,t), that scales linearly with p{ip,t), i.e. 
CipC) = Co + cip{ip,t)- We normalize this fac¬ 
tor in simh a way that its mean global annual 
value is = 1. Thanks to the normaliza¬ 

tion condition, it is possible to calculate the pa¬ 
rameters Co and ci in terms of a single parame¬ 
ter, IR = max{^((/j,t)} /min t)}, which rep¬ 
resents the ratio between the maximum and mini¬ 
mum values of C at any latitude and orbital phase 


(see Vladilo et al. 2013 §A.2.1). With the adop¬ 


tion of the modulation term, the complete expres¬ 
sion for the transport coefficient becomes 


J-^o ^dry,o 


1 -^o * (§md/Smd,o) 

1 + Ac 


• (26) 


The mean global annual value of this expression 


equals Eq. (25) thanks to the normalization con¬ 


dition C((/5,t) = 1. This formalism introduces a 
dependence on t and on the axis obliquitjj^ in the 
transport coefficient. 

Empirical support for the adoption of the 
modulation term comes from the im¬ 

proved match between the observed and predicted 
temperature-latitude profile of the Earth. In the 
left panel of Fig. [^we show that it is not possible 
to accurately match the Earth profile by varying 
Do at constant = 1 (i-®- IR = !)• This 

is because the whole profile becomes flatter with 
increasing Do and the values of Do sufficiently 
high to provide the desired smooth temperature 
distribution inside the tropics yield a profile which 
is too flat in the polar regions. This problem can 
be solved with the introduction of the modulation 
factor (. In the right panel of Fig. [^we show that 
by increasing !R the profile declines faster at the 
poles while becoming slightly flatter at the equa¬ 
tor. This behavior is different from that induced 
by changes of Do and provides an extra degree of 
freedom to match the observed profile. For the 
time being, the parameter IR can be tuned to fit 
the Earth model, but cannot be validated with 
other planets. The validation of IR in rocky plan¬ 
ets different from the Earth could be addressed by 
future GCM calculations. Meantime, the uncer¬ 
tainty related to the choice of this parameter in 
other planets can be estimated by repeating the 
climate simulations for different values of IR. Given 
the lack of solid theoretical support for the adop¬ 
tion of the formalism, it is safe to use the 

smallest possible value of IR (i.e. closest to unity) 
that allows the Earth profile to be reproduced. 
With the upgraded calibration of the Earth model 
presented here (Appendix B) we have been able 
to adopt a lower value (IR = 2.2, Table [Tj ) than in 
Paper I (IR = 6, Table 2 in Vladilo et al. (2013)). 


^'^The mean diurnal value of is a function of the axis 

obliquity. 
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2.1.3. Ocean transport 


The algorithm that describes the energy trans¬ 
port has been derived assuming that most of the 
meridional transport is performed by the atmo¬ 


sphere rather than the ocean (12.1). This is a rea¬ 


sonable assumption in the Earth climate regime, 
where the atmosphere contributes 78% of the total 
transport in the Northern Hemisphere and 92% in 
the Southern Hemisphere at the latitude of max¬ 


imum poleward transport (Trenberth & Caron 


2001). In order to assess the importance of the 


ocean contribution in different planetary regimes 
one needs to run GCM simulations featuring the 
ocean component. This is a difficult task because 
the ocean circulation is extremely dependent on 
the detailed distribution of the continents and be¬ 
cause the time scale of ocean response is much 
longer than that of the atmosphere. As a re¬ 
sult, one should run GCMs with a detailed de¬ 
scription of the geography for a large number of 
orbits in order to include the ocean transport in 
the modellization of exoplanets. With this type 
of climate simulation it would be impossible to 
perform an exploratory study of exoplanet surface 
temperature, which is the aim of our model. Not 
to mention the fact that the choice of a detailed 
description of the continental distribution in ex¬ 
oplanets is completely arbitrary. It is therefore 
desirable to find simplified algorithms able to in¬ 
clude the ocean transport in zonal models, such 
as the ESTM. To this end, one should perform 3D 
numerical experiments aimed at investigating how 
the energy transport is partitioned between the 
atmosphere and the ocean in a variety of plan¬ 
etary conditions. Preliminary work of this type 
suggests that the energy transport of wind-driven 
ocean gyrevary in a roughly similar fashion to 
the energy transport of the atmosphere as exter¬ 
nal parameters vary ( Vallis fc Farneti||2009 1. The 
existence of mechanisms of compensation that reg¬ 
ulate the relative contribution of the atmosphere 


and the ocean to the total transport ( 

Bjerknes 

1964| IShaffrey & Sutton 20061 

van der Swaluw 

et al. |2007 Lucarini & Ragone 

2011 ) may also 


help build a simplified description of the atmo¬ 
sphere/ocean transport. In the case of the Earth, 
we note that the total transport is remarkably sim- 


In oceanography the term gyre refers to major ocean cir¬ 
culation systems driven by the wind surface stress. 


see 


ilar in the Southern and Northern hemisphere 
Fig. 0 in spite of significant differences between 
the two hemispheres in terms of the relative con¬ 


tribution of the ocean and atmosphere (e.g. Tren- 
berth &: Garon|[ 2 MI Fig. 7). 


2.2. The vertical transport 


The outgoing longwave radiation and the top- 
of-atmosphere albedo are parametrized using sin¬ 
gle atmospheric column calculations. In the 
present version of the ESTM, the single column 
calculations are performed with standard radi¬ 
ation codes developed at the National Genter 
for Atmospheric Research (NGAR), as part of 
the Gommunity Glimate Model (GGM) project 
NGAR-GGM (Kiehl et al. I998| ). To access these 
codes we use the set of routines GliMT (Pierre- 


humbert]|2010 Gaballero||2012 ). 


The GGM code employs an Earth-like atmo¬ 
spheric composition, with the possibility to change 
the amount of non-condensable greenhouse gases 
(i.e. GO 2 and GH 4 ). We adopt pG 02 = 380ppmv 
and PGH 4 = 1 . 8 ppmv as the reference values for 
the Earth’s model. These values can be changed 
as long as they remain in trace abundances, as in 
the case of the Earth. The relative humidity, q, 
is fixed to limit the huge amount of calculations 
and the dimensions of the tables described below. 
We adopt q — 0.6, a value consistent with the 
global relative humidity measured on Earth. A 
low effective humidity {q ^ 0 . 6 ) is predicted self- 
consistently by 3-D dynamic climate models as a 
result of subsidence in the Hadley circulation (e.g. 
Ishiwatari et al. Adoption of saturated wa¬ 

ter vapour pressure {q = I) tends to understimate 
the OLR at high temperatures, leading to exces¬ 
sive heating of the planet. 


2.2.1. Outgoing long-wavelength radiation 

We use a column radiation model scheme for 
a cloud-free atmosphere to calculate the Outgo¬ 
ing Long-wavelength Radiation (OLR), i.e. the 
thermal infrared emission that cools the planet. 
The OLR calculations are repeated a large num¬ 
ber of times in order to cover a broad inter¬ 
val of surface temperature, T, background pres¬ 
sure, p, gravity acceleration, g, and partial pres¬ 
sure of non-condensable greenhouse gases. The 
results of these calculations are stored in tables 










































Fig. 2.— Annual temperature-latitude profile of the Earth model obtained by varying the transport pa¬ 
rameters Do and Jl. Left panel: variation of Do at IR = 1 (i.e. without seasonal modulation of the transport 


coefficient). Right panel: variation of the strength of seasonal modulation, IR, at constant Do- See 12.1.2 


OLR=OLR(T,p, 5 ,pC 02 ,pCH 4 ). In the course of 
the simulation, these tables are interpolated at 
the zonal and instantaneous value of T = T{}p,t). 
The long-wavelength forcing of the clouds is sub¬ 
tracted at this stage, taking into account the zonal 
cloud coverage, as we explain in §2.3.2[ The total 
CPU time required to cover the parameter space 
(T,p,g,pC 02 ,pCH 4 ) is relatively large. However, 
once the tables are built up, the simulations are 
extremely fast. 

2.2.2. Incoming short-wavelength radiation 

The top-of-atmosphere albedo, A, is calculated 
with the CCM code, to take into account the 
transfer of short-wavelength stellar photons in the 
planet atmosphere. In each atmospheric column 
we calculate the fraction of stellar photons that is 
reflected back in space for different values of T, p, 
g, pC 02 , surface albedo, Og, and zenith distance 
of the star, Z. In practice, for each set of values 
{g, pC 02 , PCH 4 ), we calculate the temperature 
and pressure dependence of A. Then, for each set 
of values {T,p, g,pC 02 ,pCH 4 ) the calculations are 
repeated to cover the complete intervals of surface 
albedo, 0 < Og < 1 , and zenith distance, 0 ° < 
Z < 90°. The results of these calculations are 
stored in multidimensional tables. In the course 
of the ESTM simulations these tables are inter¬ 
polated to calculate A as a function of the zonal 
and instantaneous values of (T,p, g,pC 02 , Os, Z). 


Each single column calculation of A is relatively 
fast, compared to the corresponding calculation 
of /. However, due to the necessity of covering 
a larger parameter space, the preparation of the 
tables A = A{T,p,g,pC02,asT Z) requires a com¬ 
parable CPU time. 


2.2.3. Caveats 


The CCM calculations that we 


pressure broadening (Kiehl et al. 1998 


use include 
and refs. 


therein), but not collision-induced absorption. As 
a result, the model may underestimate the atmo¬ 
spheric absorption at the highest values of pres¬ 
sure. To avoid physical conditions not considered 
in the calculations we limit the surface pressure at 
p < 10 bar. 

The calculations are valid for a solar-type spec¬ 
tral distribution. The spectral type of the central 
star affects the vertical transport because of the 
wavelength dependence of the atmospheric albedo 
(e.g., Selsis et al.|2007 ). The present version of the 
ESTM should be applied to planets orbiting stars 
with spectral distributions not very different from 
the solar one. 
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2.3. Surface and cloud properties 

2.3.1. Zonal coverage of oceans, lands, ice and 
clouds 

The zonal coverage of oceans is a free parame¬ 
ter, fo, that also determines the fraction of conti¬ 
nents, fi = 1 — fo- In this way, the planet geogra¬ 
phy is specified in a schematic way by assigning a 
set of fo values, one for each zone. The zonal cov¬ 
erage of ice and clouds is parametrized using algo¬ 
rithms calibrated with Earth experimental data. 
Following WK97, the zonal coverage of ice is a 
function of the mean diurnal temperature. 


fi{T) = maxjo. 


1 — e 


(T-273.15K)/10K 


}. (27) 


One problem with this formulation is that the ice 
melts completely and instantaneously as soon as 
T > 273.15 K. To minimize this effect, we intro¬ 
duced an algorithm that mimics the formation of 
permanent ice when a latitude zone is below the 
freezing point for more than half the orbital pe¬ 
riod. In this case, we adopt a constant ice cover¬ 
age for the full orbit, fi = fi{T), where T is the 
mean annual zonal temperature. 

As far as the clouds are concerned, we adopt 
specific values of zonal coverage for clouds over 
oceans and continents. The dependence of the 
cloud coverage on the type of underlying surface 
has long been known (e.g. Kondrate^ 19691 and 


has been quantified in recent studies (e.g. Sanroma 


& Palle 

2012 

Stubenrauch 

2013). Based on the 

results obtained by 

Sanroma & Palle (2012 

), we 


adopt 0.70 and 0.60 for the cloud coverage over 
oceans and lands, respectively. In this way, the 
reference Earth model (Appendix B) predicts a 
mean annual global cloud coverage <fc,o> = 0.67, 
in excellent agreement with most recent Earth 
data ( Stubenrauch|20i3 ). With our formalism the 
cloud coverage is automatically adjusted for plan¬ 
ets with cloud properties similar to those of the 
Earth, but different fractions of continents and 
oceans. Since the coverage of ice, fi, depends on 
the temperature, the model simulates the feedback 
between temperature and albedo. 

2.3.2. Cloud radiative properties 

The albedo and infrared absorption of the 
clouds have cooling and warming effects of the 
planet surface, respectively. Even with specifically 


designed 3D models it is hard to predict which of 
these two opposite effects dominate. The single¬ 
column radiative calculations used in studies of 
habitability usually assume cloud-free radiative 
transfer and tune the results by playing with the 


albedo ( 

Kasting 

1988 

'Casting et al. 

1993 

Kop- 

parapu et al. 

2013 

201 ^ 

;). The approach that we 


adopt with the ESTM is to parametrize the albedo 
and the long-wavelength forcing of the clouds as¬ 
suming that their global properties are similar to 
those measured in the present-day Earth. Follow¬ 
ing WK97, we express the albedo of the clouds 
as 

Oc = a + j3Z (28) 

where the parameters a and j5 are tuned to fit 
Earth experimental data of cloud albedo as a func¬ 
tion of stellar zenith distance (Cess 1976). For 


clouds over ice, we adopt the same albedo of frozen 
surfaces (see Table . To take into account the 
long wavelength forcing of the clouds, we subtract 
<OLR>ci,o (/c/</c.o>) from the clear-sky OLR 
obtained from the radiative calculations, where 
<OLR>ci,o = 26.4 W m“^ is the mean global 
long wavelength forcing of the clouds on Earth 
( Stephens et al.|2()T2 ), fc is the mean cloud cover¬ 
age in each latitude zone, and <fc,o> = 0.67 the 
mean global cloud coverage of the reference Earth 
model. 

The fact that the ESTM accounts for the mean 
radiative properties of the clouds is an improve¬ 
ment over classic EBMs, but one should be aware 
that the adopted parameterization is only valid 
for planets with global cloud properties similar to 
those of the Earth. This is a critical point because 
the cloud radiative properties may change with 
planetary conditions, as suggested by 3D simula¬ 


tions of terrestrial planets (e.g. Leconte et al. 2013 


Yang et al. 2013). To some extent, we can simu¬ 


late this situation by changing the ESTM cloud¬ 
forcing parameters. An example of this exercise is 
provided in Fig. If the predictions of 3D exper¬ 
iments become more robust, it could be possible in 
the future to introduce a new ESTM recipe for ex¬ 
pressing the cloud forcing as a function of relevant 
planetary parameters. 

2.3.3. The surface albedo 

The mean surface albedo of each latitude zone is 
calculated by averaging the albedo of each type of 
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surface present in the zone, weighted according to 
its zonal coverage. For the surface albedo of conti¬ 
nents and ice we adopt the fiducial values listed in 
Table [TJ The albedo of the oceans is calculated as 
a function of the stellar zenith distance, Z, using 
an expression calibrated with experimental data 
( Briegleb et al.|[l9^ Enomoto]|2007 1 


0.026 


tto = 


(1.10.065) 

-b 0.15(^ - 0.1) {[I - 0.5) {[I - 1.0) , (29) 

where ^ = cosZ. Also clouds are treated as 
surface features, with zonal coverage and albedo 


parametrized as explained above ((2.3.2). 


2.3.4- Thermal capacity of the surface 

The term C is calculated by averaging the ther¬ 
mal capacity per unit area of each type of surface 
present in the corresponding zone according to its 


zonal coverage ((2.3.1). The parameters used in 


these calculations are representative of the ther¬ 
mal capacities of oceans and solid surface (Table 
[^. For the reference Earth model the ocean con¬ 
tribution is calculated assuming a 50 m, wind- 
mixed ocean layeip^ (Williams & Kasting 


1997 


Pierrehumbert 2010). The atmospheric contribu¬ 


tion is calculated as 


Ca 


a 


atm,o 


(30) 


where Cp and p are the specific heat capacity 
and total pressure of the atmosphere, respectively 
( PierrehumbertpOlO I. The atmospheric term en¬ 
ters as an additive contribution to the ocean and 
solid surface terms. Its impact on these parame¬ 
ters is generally small, the ocean contribution be¬ 
ing the dominant one. 

The strong thermal inertia of the oceans im¬ 
plies that the mean zonal C has an “ocean-like” 
value even when the zonal fraction of lands is com¬ 


parable to that of the oceans (Williams & Kast- 
ing||1997 ). This weak point of the longitudinally- 
averaged model can be by-passed by adopting an 
idealized orography with continents covering all 
longitudes (see (C.6). 


The thermal capacity of the ocean can be changed to simu- 
late idealized aquaplanets with a thin layer of surface water 


(e.g., [3.11. 


2.4. The insolation term S 

The zonal, instantaneous stellar radiation S = 
S{ip,t) is calculated from the stellar luminosity, 
the keplerian orbital parameters and the inclina¬ 
tion of the planet rotation axis. The model calcu¬ 
lates S also for eccentric orbits. Details on the im¬ 


plementation of S can be found in Paper I (Vladilo 


et al. 2013, §A.5). At variance with that paper, the 
ESTM takes also into account the vertical trans¬ 


port of short-wavelength photons (see (2.2.2). 


2.5. Limitations of the model 

In spite of the above-mentioned improvements 
over classic EBMs, the adoption of the zonal en¬ 
ergy balance formalism at the core of the ESTM 
leads to well known limitations intrinsic to EBMs. 
One is that zonally averaged models cannot be ap¬ 
plied to tidally-locked planets that always expose 
the same side to their central star: such cases re¬ 


quire specifically designed models (e.g., Kite et al. 


2011 |Menou|2013[ [Mills fc Abbot|2013[[Yang et al. 


2013). Also, it should be clear that the ESTM 


does not track climatic effects that develop in the 
vertical direction, even though the atmospheric 
response is adjusted according to latitudinal and 
seasonal variations of T and Z. In spite of these 
limitations, the EBM at the core of our climate 
tools provides the flexibility that is required when 
many runs are needed or when one wants to com¬ 
pare the impact of different parameters uncon¬ 
strained by the observations. At present time this 
is still unfeasible with GCMs and even with Inter¬ 
mediate Complexity Models. While GCMs are in¬ 
valuable tools for climate change studies on Earth, 
they are heavily parameterized on current Earth 
conditions, and their use in significantly different 
conditions raises concern. In particular, the paper 


by Stevens & Bony (2013) caused serious worries 
about the use of GCMs in “unconstrained” situ¬ 
ations such as those encountered in habitability 
studies. 

3. Model calibration and validation 

The ESTM is implemented in two stages. First, 
a reference Earth model is built up by tuning the 
parameters to match the present Earth climate 
properties (see Appendix B). Then we use results 
obtained from 3D climate experiments to tune pa¬ 
rameters or validate algorithms that are meant to 
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be applied in Earth-like planets. Here we present 
a test of validation of the algorithms that describe 
the meridional transport. This test is a concrete 
example of how results obtained by GCMs can be 
used to validate the model. 


3.1. Validation of the meridional trans¬ 
port 

To perform this test we used a study of the at¬ 
mospheric dynamics of terrestrial exoplanets per¬ 


formed by Kaspi & Showman (2014). These au¬ 


thors employed a moist atmospheric general cir¬ 
culation model to test the response of the at¬ 
mospheric dynamics over a wide range of planet 
parameters. Specifically, they used an idealized 
aquaplanet with surface covered by a uniform 
slab of water 1 m thick; only vapor-liquid phase 
change was considered; the albedo was fixed at 
0.35 and insolation was imposed equally between 
hemispheres; the remaining parameters were set 
to mimic an Earth-like climate. To validate the 
ESTM with the results found by [Kaspi &: Show-] 


man 


(2014) we modified the Earth reference model 
as follows. The axis obliquity was set to zero; 
the temperature-ice feedback was excluded; the 
albedo was fixed at A = 0.35; the fraction of 
oceans was set to 1 , adopting a thermal capac¬ 
ity corresponding to a mixing layer 1 m thick. 
With this idealized planet model we performed 
several sets of simulations, varying the planet ro¬ 
tation rate, surface flux, radius, and surface pres¬ 
sure. To validate the ESTM we analyze the 
mean annual equator-to-pole temperature differ¬ 
ence, ATep, which is critical for a correct esti¬ 
mate of the latitude temperature profile and of 
the surface habitability. The results of the tests 
are shown in Fig. where we compare the ATpp 
values predicted by the 3D model (diamonds) with 
those obtained from the ESTM (solid lines). We 
also plot the predictions of a “dry” transport 
model (dashed lines) obtained by setting A = 0 
in Eq. ( |20| . Finally, for the sake of comparison 
with previous EBMs, we plot the results obtained 
from a “basic ESTM” without moist term (A = 0) 
and without diabatic forcing ternf^ (dotted lines). 
In using this “basic” model, we test some alterna¬ 
tive scaling laws for the parameterization of the 


rotation rate, surface pressure and radius, as we 
explain below. 


3.1.1. Rotation rate 


In this experiment all parameters were fixed, 
with the exception of the planet rotation rate, 
H, that was gradually increased from 1/10 to 10 
times the Earth value, H 0 . The results of this 
test are shown in the top-left panel of Fig. One 
can see that the ESTM and GCM results show 
a similar trend, with a good quantitative agree¬ 
ment at O > 0.3 00 , but not at low rotation rate. 
This result is expected since our parameterization 
is appropriate to simulate planets with horizon¬ 
tal transport dominated by mid-latitude eddies, 
i.e. planets with relatively high rotation rate (see 


(2.1 1 . The dotted line in this figure shows the re¬ 


sults of the basic model obtained by replacing the 
term 0“^^® in Eq. (23) with the strong er depen- 
dence 0“^ adopted in previous work (e.g. Williams 


& Kasting 1997 Vladilo et al. 2013| ). One can 
see that this strong dependence on rotation rate 
is not supported by the 3D model, while the more 
moderate dependence D oc adopted in the 

ESTM yields a much better agreement with the 
GCM experiments. 


3.1.2. Stellar flux 

In the top-right panel of Fig. [^we show the re¬ 
sults obtained by varying the insolation from 100 
to 2000 Wm“^, i.e. from 0.07 to 1.47 times the 
present-day Earth’s insolation. The behavior pre¬ 
dicted by the 3D model is bimodal, with a rise of 
ATep up to an insolation of ~ 800 Wm“^ and a 
decline at higher values of stellar flux. According 
to Kaspi & Showman (2014) the decline is trig¬ 


gered by the rise of the moist transport efficiency 
resulting from the increase of temperature and wa¬ 
ter vapor content. The moist ESTM is able to cap¬ 
ture this bimodal behavior, even though a reason¬ 
able agreement with the 3D experiments is only 
found in a range of insolation ±25% around the 
present-Earth’s value (shaded area in the figure). 
The dry model (dashed line) is unable to capture 
the bimodal behavior of ATep versus flux. The 
basic model is even more discrepant (dotted line). 


To ignore the diabatic forcing term we set {ASR} = 1 
in the scaling law ||23|l. 
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Fig. 3.— Predicted values of equator-to-pole temperature difference, ATep, versus planet rotation rate 
(top left), stellar flux (top right), surface atmospheric pressure (bottom left), and planet radius (bottom 
right) obtained from different climate simulations of an Earth-like aquaplanet with axis obliquity e = 0°, 
fixed albedo {A = 0.35), and no ice. Magenta diamonds: 3D GCM simulations (Kaspi & Showman 2014| ). 
Solid and dashed curves: ESTM (moist and dry transport, respectively). Dotted curves: ESTM dry model 
without diabatic forcing term and with different prescriptions for D (see legend in each panel and (3.1 ) The 
shaded area in the top right panel brackets a range of ±25% the present-day value of Earth’s insolation. 
The shaded area in the bottom right panel brackets the range of masses 0.1 to 8 M® for rocky planets with 
mean density p = p®. 
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3.1.3. Surface pressure or atmospheric columnar 
mass 


In the bottom-left panel of Fig. we show 
the results obtained by varying the surface pres¬ 
sure p of the idealized aquaplane! from 0.2 to 
20 bar. Since the surface gravity is not varied, 
this experiment is equivalent to vary the atmo¬ 
spheric columnar masj^ p/g, from 0.2 to 20 times 
that of the Earth. Theoretical considerations indi¬ 
cate that the efficiency of the horizontal transport 
must increase with increasing p/g [e.g. Eq. ([^], 
and equator-pole temperature differences should 
decrease as a result. The 3D model predicts a 
monotonic decrease of ATep, in line with this ex¬ 
pectation. However, the decrease is milder than 
expected by the basic model with a simple law 
D (X p/g (dotted line). The models with diabatic 
forcing (solid and dashed lines) predict a more 
moderate decrease, D oc {p/gY^^ [Eq. (23)], and 
are in much better agreement with the results of 
the 3D experiments. The agreement of the moist 
ESTM (solid line) is remarkable in the range of 
high columnar mass. 


3 . 1 . 4 . Planet radius or mass 

In this experiment all planetary parameters, in¬ 
cluding the columnar mass p/g, are fixed while 
changing the planet radius. Assuming a constant 
mean density, p = p^, this is equivalent to scale 
the planet mass as M oc R^. The results are 
shown in the bottom-right panel of Fig. where 
3D models predict an increase of ATpp with in¬ 
creasing radius and mass, indicating that the hor¬ 
izontal transport becomes less efficient in larger 
planets. This is in line with theoretical expecta¬ 
tions which suggest that the transport coefficient 
decreases with increasing radius, possibly with a 
quadratic law [e.g. Eq. ©I- However, the in¬ 
crease of ATep appears to be too sharp if we adopt 
the basic model with D oc R~^ (dotted line). The 
models with diabatic forcing (solid and dashed 
lines) predict a moderate decrease, D oc 
[Eq. (13], in line with the 3D predictions. In 


^ In this experiment, Kaspi Showman i 20141 adopted a 
constant optical depth of the atmosphere to focus on hor¬ 
izontal transport, rather than vertical transfer effects. For 
the sake of comparison with their experiment, we used a 
constant value of atmospheric columnar mass in the OLR 
and TOA-albedo calculations, while changing p/g in the 
diffusion term. 


the range of masses typical of terrestrial planets 
(shaded area in the bottom-right panel of Fig. 
the predictions of the ESTM are very similar to 
those obtained by Kaspi & Showman (2014). 


4. Applications 

After the calibration and validation, we apply 
the model to explore the dependence oiT{ip, t) and 
the mean global surface temperaturf^ T, on a 
variety of planet parameters. At variance with the 
validation tests, we now consider all the features 
of the model, including the ice-albedo feedback.In 
Appendix C we present simulations of idealized 
Earth-like planets. Here we describe a test study 
of exoplanet habitability. 


4.1. Exoplanets 

The modelization of the surface temperature of 
exoplanets is severely constrained by the limited 
amount of observational data. Typically, one can 
measure the stellar and orbital parameters and a 
few planetary quantities, such as the radius and/or 
mass. From the stellar and orbital data one can 
estimate the planet insolation and its seasonal evo¬ 
lution. From the radius and mass one can estimate 
the surface gravity which enters in the parameter¬ 
ization of the atmospheric columnar mass. Unfor¬ 
tunately, many planet quantities that are required 
for the modelization are currently not observable. 
These include the atmospheric compositioiFI sur¬ 
face pressure, ocean/land distribution, axis obliq¬ 
uity and rotation period. Taking advantage of the 
flexibility of the ESTM, we can perform a fast ex¬ 
ploration of the space of the unknown quantities, 
treating them as free parameters. From the ap¬ 
plication of this methodology we can assess the 
relative importance in terms of climate impact of 
the planet quantities that are not measurable. In 
addition, we can constrain the ranges of parame¬ 
ters values that yield habitable solutions. 

We show two examples of application of this 
methodology. First we consider a specific exo¬ 
planet chosen as a test case, then we introduce 
a statistical ranking of planetary habitability. We 


^®The average in latitude is weighted in area and the average 
in time is performed over one orbital period. 

^®At present time it is possible to measure the atmospheric 
composition of selected giant planets, but not of Earth-size 
terrestrial planets. 
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Kepler62e Kepler62e 




Fig. 4.— Influence of planetary parameters not constrained by observations on the mean global surface 
temperature, T, of Kepler62-e. Symbols: predictions of ESTM simulations plotted as a function of surface 
pressure, p, for all possible combinations of rotation rate, axis obliquity and ocean fraction described in 
§4.1.1[ Circles: habitable solutions. Diamonds: solutions with water vapor column above the critical limit 
discussed in SQ The symbol size scales with the fractional habitability, /iiw Left: Earth atmospheric 
composition. Right: Earth-like atmospheric composition with a one hundred-fold increase of pC02. Dashed 
line: equilibrium temperature of Kepler62-e ( ]Borucki]|2013 ). Dotted line: water boiling point as a function 
of p. 
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adopt an index of habitability, hiw, based on the 
liquid water criteriorp^ 


4-. 1.1. Kepler-62e as a test case 


The test-case exoplanet was chosen using three 
criteria. The first is that the planet should be 
of terrestrial type, i.e. rocky and without an ex¬ 
tended atmosphere, in order to be suitable for 
the application of the ESTM. We used the radius 
for a preliminary characterization of the planet, 
since evidence is accumulating for the existence of 
a gradual transition, correlated with radius, be¬ 
tween planets of terrestrial type and planets with 


rocky cores but extended gas envelopes (e.g. Wu 
fc Lithwi'^|2013 Marcy|[2M4 1 . We restricted our 
search to planets with R < l.Ti?®, the thresh¬ 
old for terrestrial planets found in a statistical 
study of size, host-star metallicity and orbital pe¬ 
riod ( Buchhave et al.|2014 ). As a second criterion, 
we required the orbital semimajor axis to be larger 
than the tidal lock radius, since the ESTM cannot 
be applied to tidally locked planets. Einally, given 
the extreme dependence of habitability on insola¬ 


tion (see e.g. Fig. 10 1 , we selected planets with an 
insolation within ±50% of the present-day Earth 
value. By querying the Exoplanet Orbit Database 


(Wright 20111 at exoplanets.org, we found that 
only Kepler-62 e ( Borucki|2013 1 satisfies the above 
criteria. The radius, R = l.Gli?®, and orbital 
period, P = 122 d, suggest that Kepler-62 e is 


probably of terrestrial type (see Buchhave et al. 


2014, Fig. 2). Its insolation is only 19% higher 


than the Earth’s value, and its semimajor axis, 
a = 0.427 AU, is larger than the tidal lock ra- 
diuj^ rti = 0.31 AU. 

To run the ESTM simulations of Kepler-62e we 
adopted at face value the radius, semimajor axis, 
eccentricity and stellar flux provided by the ob¬ 


servations (Borucki 2013). Unfortunately, only a 


loose upper limit (M < 36 M 0 ) is available for 


We define a function = 1 when T{ip, t) is inside the 

liquid-water temperature range; = 0 outside the 

range ( [Spiegel et al.|2008| jVladilo et a.l.|2013| l . The index 
is the global and orbital average value of 
The tidal lock radius was calculated with the expression 
rxL = 0.027{Pot/Q)^l^Ml^^ (Kasting 1993), where Po is 
the original rotation period, t is the amount of time during 
which the planet has been slowed down, Q~^ is the specific 
dissipation function, and M* is the stellar mass; we adopted 
Po = 0.5 days, t = 10® yr, and Q = 100. 


the mass, so that the surface gravity g is poorly 
constrained at the present time. For illustrative 
purposes, we adopt g = 1.5 (M = 3.9 M 0 ), cor¬ 
responding to a mean density 5.1gcm“^, similar 
to that of the Earth = 5.5gcm“^). As far as 
the atmosphere is concerned, we vary the surface 
pressure in the range p G (0.03, 8 ) bar and the CO 2 
partial pressure in the range pCO 2 /(pCO 2)0 G 
(1,100). We adopt 3 representative values of ro¬ 
tation rate, n/ll 0 G (0.5,1,2), axis obliquity, 
e G (0°, 22.5°, 45°), and ocean coverage, fo G 
(0.5,0.75,1.0). For the remaining parameters we 
adopt the Earth’s reference values. For each value 
of CO 2 partial pressure we run simulations cover¬ 
ing all possible combinations of background pres¬ 
sure, rotation rate, axis obliquity and ocean cov¬ 
erage listed above. Part of the results of these 
simulations are shown in Figs. |^and[^ 

^ In Fig. I^we plot the mean global temperature, 
T, obtained for two different values ofpC 02 , spec¬ 
ified in the legends. At each value of p, we show 
the values of T obtained from all possible combi¬ 
nations of n, e, and fo. The typical scatter of T 
due a random combination of these 3 parameters 
is ~ 10-20 K. On top of this scatter, the most re¬ 
markable feature is a positive trend of T versus p 
extending over an interval of « 60 K. Within the 
limits of application of the model, these results in¬ 
dicate that an uncertainty on p within the range 
expected^Jor terrestrial planet^^ has stronger ef¬ 
fects on T than uncertainties of rotation rate, axis 
obliquity and ocean coverage. In fact, variations 
of p have strong effects both on T{tf, t) and T be¬ 
cause they are equivalent to variations of atmo¬ 
spheric columnar mass, p/g, which affect both the 
latitudinal transport (i.e. the surface temperature 
distribution), and the radiative transfer (i.e. the 
global energy budget). 

The results shown in Fig. constrain the in¬ 
terval of p that allows Kepler-62e to be habitable. 
At Jiigh p, the habitability is limited by the rise 
of T, which eventually leads to a runaway green¬ 
house instability. The red diamonds in Fig. 
indicate cases where the water vapor column ex¬ 
ceeds the critical value that we tentatively adopt 
as a limit for the onset of such instability ()fA|. 
At low p, two factors combine to limit the hab- 


^® Mars and Venus have a surface pressure of ~ 0.006 bar and 
~ 90 bar, respectively. 
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itability. One is the onset of large temperature 
excursions and the other the decrease of the wa¬ 
ter boiling point (dotted lines in Fig. |^. As a 
result, the fraction of planet surface outside the 
liquid water temperature becomes larger at low p. 
To evidentiate this effect, we have scaled the size 
of the symbols in Fig. according to the value of 
h\^. One can see that h\^ tends to become smaller 
at low p, especially when T approches the tem¬ 
perature regime where the ice-albedo feedback be¬ 
comes important. In some cases, not shown in the 
figure, the planet undergoes a complete snowball 
transition and the habitability becomes zero. The 
effect of temperature excursions is shown in Fig. 
where we plot as a function of p the average value 
of ATep obtained from all possible combinations 
of f2, e, and fo- One can see that the temperature 
excursions become large with decreasing level of 
CO 2 ; this happens because at low CO 2 the tem¬ 
perature is sufficiently low for the development of 
the ice-albedo feedback, and because the lowered 
IR optical depth of the atmosphere is less effective 
in reducing the effect of the geometrically-induced 
meridional insolation variation at the surface. 


Fig. a shows that the equilibrium temperature 
of Kepler-62e, = 270 ± 15K ( |Borucki|[M^ 
dashed line), lies at the lower end of the predicted 
T values. The difference T — Teq increases with 
atmospheric columnar mass because the estimate 
of Teq does not consider the greenhouse effect. 


4 . 1 . 2 . Statistical ranking of planetary habitability 


By performing a large number of simulations for 
a wide combination of parameters we can quantify 
the habitability in a statistical way. To illustrate 
this possibility with an example we consider again 
the test case of Kepler-62e. We perform a statis¬ 
tical analysis of the results obtained from all com¬ 
binations of n, e, and fo values adopted at a given 
p. We tag as “non habitable” the cases with a 
snowball transition and those with a critical value 
of water vapor. For each set of parameters {17, e, 
fo} we count the number of cases that are found 
to be habitable, n^, over the total number of sim¬ 
ulations, Uf. From this we calculate the fraction 
iph = nfi/nt, which represents the probability for 
the planet to be habitable if the adopted parame¬ 
ter values are equally plausible a priori. We then 
call = (1/n/i) X)r=i(^iw)i the mean surface 



Fig. 5.— Average equator-pole temperature dif¬ 
ference, <ArEp>, obtained from ESTM simula¬ 
tions of Kepler62-e (j |4.1.1 1 , plotted as a function 
of surface pressure, p. Each curve is calculated at 
constant PCO 2 , as specihed in the legend. 



Fig. 6.— Ranking index of habitability, r?i, ob¬ 
tained from ESTM simulations of Kepler-62e and 
an Earth twin (1 ]4.1.2 ) plotted as a function of sur¬ 
face pressure, p. Symbols for Kepler-62e as in Fig. 
m Crossed circles: Earth twin. 
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habitability of the nh sets that yield a habitable 
solution. At this point we define a “ranking index 
of habitability” 

^ nh 

'^h — '4’h ^ ^ ' (hlw)2 ■ (^1) 

nt ^ 

This index takes into account at the same time the 
probability that the planet is habitable and the av¬ 
erage fraction of habitable surface. As an example, 
in Fig. [^we plot Vh versus p for the different values 
of pC02 considered in our simulations of Kepler- 
62e. One can see that rh — ^ only in a limited 
range of surface pressure. At very low pressure, 
the index becomes lower because the fraction of 
habitable surface decreases and because drops 
below unity when a snowball transition is encoun¬ 
tered. At high pressure iph drops when the water 
vapor limit is encountered. As an example of ap¬ 
plication, we can constrain the interval of p suit¬ 
able for the habitability of Kepler62-e. For pC02= 
(PCO 2 )© (triangles) the requirement of habitabil¬ 
ity yields the limit p ^ 3 bar. As PCO 2 increases 
(squares and pentagons), the upper limit becomes 
more stringent {p < Ibar), but the planet has a 
higher probability of being habitable at relatively 
low pressure. 

Clearly, the index r/j does not have an absolute 
meaning since its value depends on the choice of 
the set of parameters. However, by choosing a 
common set, the index rh can be used to rank the 
relative habitability of different planets. As an 
example, in Fig. [^we plot versus p obtained for 
an Earth twirp^ with the same sets of parameters 
{H, e, /o} adopted for Kepler-62e. One can see 
that at p > 1 bar the Earth twin (crossed circles) 
is more habitable than Kepler-62e for the adopted 
set of parameters, while at p < 1 bar Kepler-62 is 
more habitable than the Earth twin. 

5. Summary and conclusions 

We have assembled the ESTM set of climate 
tools (Figure 1) to model the latitudinal and sea¬ 
sonal variation of the surface temperature, T{ip, t), 
on Earth-like planets. The motivation for build¬ 
ing the ESTM is twofold. From the general point 

Here an Earth twin is a planet with properties equal to 
those of the Earth (including insolation, radius, gravity and 
atmospheric composition), but with unknown values of p, 
fl, e and fo, which are treated as free parameters. 


of view of exoplanet research, Earth-size planets 
are expected to be rather common, but difficult 
to characterize with experimental methods. From 
the astrobiological point of view, Earth-like plan¬ 
ets are excellent candidates in the quest for hab¬ 
itable environments. A fast simulation of T(ip,t) 
enables us to characterize the surface properties 
of these planets by sampling the wide parameter 
space representative of the physical quantities not 
measured by observations. The detailed modeliza- 
tion of the surface temperature is essential to es¬ 
timate the habitability of these planets using the 
liquid water criterion or a proper set of thermal 
limits of life (e.g. ClarkepOld ). 

The ESTM consists of an upgraded type of 
EBM featuring a multi-parameter description of 
the physical quantities that dominate the verti¬ 
cal and horizontal energy transport (Fig.[^. The 
functional dependence of the physical quantities is 
derived using single-column atmospheric calcula¬ 
tions and algorithms tested with 3D climate ex¬ 
periments. Special attention has been dedicated 
to improve ((2.1 1 and validate ((3.11 the descrip¬ 
tion of the meridional transport, a weak point 
of classic EBMs. The functional dependence of 
the meridional transport on atmospheric colum¬ 
nar mass and rotation rate is significantly milder 
[see Eq. (23)] than the one adopted in previous 
EBMs. The reference Earth model obtained from 
the calibration process is able to reproduce with 
accuracy the average surface temperature proper¬ 
ties of the Earth and to capture the main features 
of the Earth albedo and meridional energy trans¬ 
port (Figs. [^and|^. Once calibrated, the ESTM 
is able to reproduce the mean equator-pole tem¬ 
perature difference, ATep, predicted by 3D aqua- 
planet models (Fig. [^. 

The ESTM simulations provide a fast “snap¬ 
shot” of T{(p,t) and temperature-based indices of 
habitability for any set of input parameters that 
yields a stationary solution. The planet parame¬ 
ters that can be changed include radius, i?, sur¬ 
face pressure, p, gravitational acceleration, g, ro¬ 
tation rate, D, axis tilt, e, ocean/land coverage 
and partial pressure of non-condensable green¬ 
house gases. The approximate limits of validity 
of the present version of the ESTM can be sum¬ 
marized as follows: 0.5 < R/R® ^ 2, p < 10 bar, 
0.5 < D/D© ^ 5, e < 45°; pC 02 and pCH4 can be 
changed, but should remain in trace abundances 
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with respect to an Earth-like atmospheric compo¬ 
sition. The requirement of a relatively high rota¬ 
tion rate is inherent to the simplified treatment of 
the horizontal transport. However, the ESTM can 
be applied to explore the early habitability of slow 
rotating planets that had an initial fast rotation. 

We have performed ESTM simulations of ideal¬ 
ized Earth-like planets to evaluate the impact on 
the planet temperature and habitability resulting 
from variations of rotation rate, insolation, at¬ 
mospheric columnar mass, radius, axis obliquity, 
ocean/land distribution, and long-wavelength 
cloud forcing (Figs. [|-[T§. Most of these quan¬ 
tities can easily induce ~ 30-40% changes of the 
mean annual habitability for parameter variations 
well within the range expected for terrestrial plan¬ 
ets. Variations of insolations within ±10% of the 
Earth value can impact the habitability up to 
100%. The land/ocean distribution mainly affects 
the seasonal habitability, rather than its mean an¬ 
nual value. The impact of rotation rate is weaker 
than predicted by classic EBMs, without evidence 
of a snowball transition at H/r 20 > 3. A general 
result of these numerical experiments is that the 
ice-albedo feedback amplifies changes of T{(p,t) 
resulting from variations of planet parameters. 


We have tested the capability of the ESTM 
to explore the habitability of a specific exoplanet 
in presence of a limited amount of observational 
data. For the exoplanet chosen for this test, 
Kepler-62e, we used the stellar flux, orbital pa¬ 
rameters and planet radius provided by the ob¬ 
servations ( Borucki|2013 ), together with a surface 
gravity g = 1.5 30 adopted for illustrative pur¬ 
poses. We treated the surface pressure, _pC 02 , ro¬ 
tation rate, axis tilt and ocean coverage as free pa¬ 
rameters. We find that T increases from « 280 K 
to « 340 K when the surface pressure increases 
between p ~ 0.03 bar and 3 bar; this trend domi¬ 
nates the scatter of ~ 10-20 K resulting from vari¬ 
ations of rotation rate, axis tilt and ocean cov¬ 
erage at each value of p. We also find that the 
surface pressure of Kepler-62e should lie above 
p « 0.05 bar to avoid the presence of a significant 
ice cover and below « 2 bar to avoid the onset of a 
runaway greenhouse instability; this upper limit is 
confirmed for different values of pC 02 and surface 
gravity. These results demonstrate the ESTM ca¬ 
pability to evaluate the climate impact of unknown 
planet quantities and to constrain the range of val¬ 


ues that yield habitable solutions. The test case 
of Kepler-62e also shows that the equilibrium tem¬ 
perature commonly published in exoplanet studies 
represents a sort of lower limit to the mean global 
temperature of more realistic models. 


We have shown that the flexibility of the ESTM 
makes it possible to quantify the habitability in a 
statistical way. As an example, we have intro¬ 
duced a ranking index of habitability, r^,, that 
can be used to compare the overall habitability 
of different planets for a given set of reference pa¬ 
rameters (14.1.2). For instance, we find that at 
p < 1 bar Kepler-62e is more habitable than an 
Earth twin for the combination of rotation rates, 
axis tilts and ocean fractions considered in our 
test, whereas the comparison favours the Earth 
twin at p > 1 bar. The index Vh can be applied 
to select the best potential cases of habitable exo¬ 
planets for follow-up searches of biomarkers. 


The results of this work indicate the level of 
accuracy required to estimate the surface habit¬ 
ability of terrestrial planets. The quality of exo¬ 
planet orbital data and host star fluxes should be 
improved to measure the insolation with an accu¬ 
racy of « 1%. In spite of the difficulty of charac¬ 
terizing terrestrial atmospheres 


le.g., 


Misra et al. 


2013), an effort should be made to constrain the 


atmospheric pressure, possibly within a factor of 
two. 


We thank Yohai Kaspi for providing results 
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suggestions received from an anonymous referee 
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this work. This research has made use of the Ex¬ 
oplanet Orbit Database and the Exoplanet Data 
Explorer at exoplanets.org. We thank Rodrigo Ca¬ 
ballero and Raymond Pierrehumbert for sugges¬ 
tions concerning the use of their climate utilities. 
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A. Running the simulations 


The ESTM simulations consist in a search for a stationary solution T((/ 3 , t) of Eq. Q. We solve the spatial 
derivates with the Euler method, with the boundary condition that the flux of horizontal heat into the pole 
vanishes. The temporal derivates are solved with the Runge-Kutta method. The solution is searched for 
by iterations, starting with an assigned initial value of temperature equal in each zone, T{ip^to) = Tstart- 
Every 10 orbits we calculate the mean global orbital temperature, T. The simulation is stopped when T 
converges within a prefixed accuracy. In practice, we calculate the increment 5T every 10 orbits and stop 
the simulation when |(5T| < 0.01 K. In most cases the convergence is achieved in fewer than 100 orbits. 
After checking that the simulations converge to the same solution starting from widely different values of 
Tstart > 273 K, we adopted Tgtart = 275 K. The choice of this “cold start” allows us to study atmospheres 
with very low pressure, where the boiling point is just a few kelvins above the freezing point; in these cases, 
the adoption of a higher Tgtart would force most of the planet surface to evaporate at the very start of the 
simulation. The adoption of a lower Tstart, on the other hand, would trigger artificial episodes of glaciation. 


In addition to the regular exit based on the convergence criterion, we interrupt the simulation in presence 
of water vapor effects that may lead to a condition of non-habitability. Specifically, two critical conditions 
are monitored in the course of the simulation. The first takes place when T{ip,t) exceeds the water boiling 
point, Tb. The second, when the columnar mass of water vapouip^ exceeds 1/10 of the total atmospheric 
columnar mass (see next paragraph). In the first case, the long term habitability might be compromised due 
to evaporation of the surface water. The second condition might lead to the onset of a runaway greenhouse 
instability ( Hart|19^ Kasting|1988 ) with a complete loss of water from the planet surface. The ESTM does 
not track variations of relative humidity and is not suited to describe these two cases. By interrupting the 
simulation when one of these two conditions is met, we limit the range of application of the simulations to 
cases with a modest content of water vapor that can be safely treated by the model. 

The limit of water vapor columnar mass that we adopt is inspired by the results of a study of the Earth 
climate variation induced by a rise of insolation ( Leconte et al.|20T3 ). When the insolation attains 1.1 times 
the present Earth value, the 3D moist model of Leconte et al. (2013) predicts an energy unbalance that 
would lead to a runaway greenhouse instability. The mixing ratio of water vapor over moist air predicted 
at the critical value of insolation is ~ 0.1 (Leconte et al. Eig. 3b). The limit of water columnar mass 

that we adopt is based on this mixing ratio. 


Eor the simulations presented in this work, we adopt a grid of iV = 54 latitude zones. The orbital period is 
sampled at Ng = 48 instants of time to investigate the seasonal evolution of the surface quantities of interest 
(e.g. temperature, albedo, ice coverage). With this set-up, the simulation of the Earth model presented below 
attains a stationary solution after 50 orbital periods, with a CPU time of ^ 70s on a 2.3 GHz processor. 
This extremely high computational efficiency is the key to perform the large number of simulations required 
to cover the broad parameter space of exoplanets. 


B. The reference Earth model 


For the calibration of the reference Earth model we adopted orbital parameters, axis tilt, and rotation 
period from Cox (2000). For the solar constant we adopted So = 1360 W m“^ and g = 9.8 ms“^ for the 


surface gravity acceleration. The zonal coverage of continents and oceans was taken from Table III in WK97. 


We adopted a relative humidity q = 0.6 (see 12.2) and volumetric mixing ratios of CO 2 and CH 4 of 380 ppmV 


and 1.8 ppmV, respectively. The surface pressure of dry air, p^ry = 1-0031 x 10® Pa, was tuned to match 
the moist surface pressure of the Earth, ptot = 1.0132 x 10® Pa. The remaining parameters of the model 
are shown in Table Some parameters were fine tuned to match the mean annual global quantities of the 
northern hemisphere of the Earth, as specified in the Table. We avoided using the southern hemisphere as a 


The water columnar mass is ~ (fJ-w/n){qPvi'^)/9)t where fiw and /r are the molecular weights of water and air; q is the 
relative humidity and p%{T) the saturated partial pressure of water vapor (e.g.,[^errehumbert|20fo||. 
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Fig. 7.— Mean annual latitude profiles of surface temperature, albedo and meridional energy flux predicted 
by the reference Earth model (solid line). Top panel: the temperature profile is compared with ERA Interim 
2m temperatures (Dee 2011) averaged in the period 2001-2013 (crosses). Middle panel: the albedo profile 


is compared with CERES short-wavelength albedo (Loeb et al. 2005 2007) averaged in the same period 
(crosses). Bottom panel: the meridional energy flux profile is compared with the total (dashed line) and 
atmospheric (dotted line) profiles obtained from the EC-Earth model (Hazeleger et al.||2010|. 


reference since its climate is strongly affected by the altitude of Antartica, while orography is not included 
in the model. In Table column 3, we show the experimental data of the northern hemisphere used as a 
guideline to tune the model parameters. In column 4 of the same table we show the corresponding predictions 
of the Earth reference model. 


In Fig. [^we show the mean annual latitude profiles of surface temperature, top-of-atmosphere albedo and 
meridional energy flux predicted by the reference model (solid line). The temperature prohle is compared 
with ERA Interim 2m temperatures ( Dee||20lT ) averaged in the period 2001-2013 (crosses in the top panel). 
Area-weighted temperature differences between observed and predicted profile have an rms of 1.1 K in the 
northern hemisphere. The albedo profiles are compared with the CERES short-wavelength albedo (|Loeb 


et al.|2005 2007) averaged in the same period (crosses in the middle panel). The ESTM is able to reproduce 


reasonably well the rise of albedo with increasing latitude. This rise is due to two factors: the dependence 
of the atmosphere, ocean, and cloud albedo on zenith distance (§ p.2.2|2.3.3 ) and the increasing coverage 
of ice at low temperature ((2.3.1). The meridional flux in the bottom panel is compared with the total 
flux (dashed line) and the atmospheric flux (dotted line) obtained from EC-Earth model (]Hazeleger et al. 
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Fig. 8.— Seasonal and latitudinal variations of surface temperature (top panels) and top-of-atmosphere 
albedo (bottom panels) of the Earth. Left panels: average ERA Interim 2m temperatures ( Dee||20lT ) and 
CERES short-wavelength albedo ( Loeb et al.||2005 20071 in the years 2001-2013. Right panels: predictions 
of the reference Earth model. The zero point of the time scale is set at the spring equinox. Blank areas in 
the albedo panels indicate regions of the seasonal-latitudinal space without incident stellar flux. 


2010). In spite of the simplicity of the transport formalism intrinsic to Eq. ([^, the model is able to capture 


remarkably well the latitude dependence of the meridional transport. 


The seasonal variations of the temperature and albedo latitudinal profiles are compared with the exper¬ 
imental data in Fig. One can see that the reference model is able to capture the general patterns of 
seasonal evolution. 


Even if the reference model has been tuned using northern hemisphere data, the predictions shown in Figs, 
and 1^ are in general agreement with the data also for most the southern hemisphere, with the exception of 
Antartica. It is remarkable that the atmospheric transport in both hemispheres is reproduced well, in spite 
of significant differences in the ocean contribution between the two hemispheres (see (2.1.31. 


Once the reference model is calibrated, some of the parameters that have been tuned to fit the present-day 
Earth’s climate can be changed for specific applications of the ESTM to exoplanets. As an example, even if 
we adopt = 0.18 for the surface albedo of continents in the reference model, we may adopt lower values, 
typical of forests, or higher values, typical of sandy deserts, for specific applications. More information on 
the parameters that can be changed is given in Table 
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C. Model simulations of idealized Earth-like planets 

In this set of experiments we study the effects of varying a single planet quantity while assigning Earth’s 
values to all the remaining parameters. We consider variations of rotation period, insolation, atmospheric 
columnar mass, radius, obliquity, land distribution, and long-wavelength cloud forcing. 


C.l. Rotation rate 


In Fig. [^we show how is affected by variations of rotation rate, the left and right panel corre¬ 

sponding to the cases = 0.5110 and = 4110, respectively. One can see that the change of the surface 
temperature distribution is quite dramatic in spite of the modest dependence of the transport coefficient on 
rotation rate that we adopt, D oc [see Eq.(21)]. The mean global habitability changes from hiw = 0.94 

in the slow-rotating case to h\^ = 0.71 in the fast-rotating case. The corresponding change of mean global 
temperature is relatively small, from T = 284 K to 290 K. These results show the importance of estimating 
rather than T, in order to quantify the habitability. 

The behavior of the mean equator-pole difference, ATePi is useful to interpret the results of this test. We 
find ATep = 28K in the slow-rotating case and ATep = 80K in the fast-rotating case. This variation of 
ATep is much higher than that found for the same change of rotation rate in the case of the KS14 aquaplanet 
(top-left panel in Fig. |^. We interpret this strong variation of ATep in terms of the ice-albedo feedback, 
which is positive and tends to amplify variations of the surface temperature. This feedback is accounted for 
in the present experiment, but not in the case of the aquaplanet. These results illustrate the importance 
of using climate models with latitude temperature distribution and ice-albedo feedback in order to estimate 
the fraction of habitable surface. 

The analysis of the ice cover evidentiates another important difference between the ESTM and classic 
EBMs. With the ESTM the ice cover increases from ~ 3% at fl = 0.5Il 0 to ~ 23% at fl = 4120. This 
increase is less dramatic than the transition to a complete “snowball” state (i.e. ice cover ~ 100 %) found 
with classic EBMs at = 3 Il 0 (e.g. Spiegel et al. 2008). This difference is due to two factors. One is 
the strong dependence of the transport on rotation rate adopted in most EBMs {D oc which is not 

supported by the validation tests discussed above (top-left panel of Fig. |^. Another factor is the algorithm 
adopted for the albedo, which in classic EBMs is a simple analytical function A = A{T), while in the ESTM 
is a multi-parameter function A = A(T,p, g,pC02,as, Z) that takes into account the vertical transport of 
stellar radiation ((2.2.2). These results show the importance of adopting algorithms calibrated with 3D 
experiments and atmospheric column calculations. 


C.2. Insolation 


In Fig. 10 we show how is affected by variations of stellar insolation, the left and right panel 

corresponding to an insolation of 0.9 and 1.1 times the present-day Earth’s insolation (S'/d ~ 341 Wm“^), 
respectively. In the first case the ESTM finds a complete snowball with T = 223 K and h\^ = 0, while in 
the latter T = 306 K and hi^ = 1 with no ice cover. These results demonstrate the extreme sensitivity of 
the surface temperature to variations of insolation and the need of incorporating feedbacks in global climate 
models in order to define the limits of insolation of a habitable planet. 


Our model is able to capture the ice-albedo feedback, but is not suited for treating hot atmospheres 
and the runaway greenhouse instability. To test the limits of the ESTM we have gradually increased the 
insolation and compared our results with those obtained by the 3D model of Leconte et al.| (20131. We found 
that the ESTM tracks the rise of T with insolation predicted by the 3D model up to 3/4 ~ 365 Wm“^. At 
higher insolation, the 3D model predicts a faster rise of T due to an increase of the radiative cloud forcing. 


By decreasing the insolation with respect to the Earth’s value, the ESTM finds solutions characterized 
by increasing ice cover. At 3/4 ~ 310 Wm“^ the simulation displays a runaway ice-albedo feedback that 
leads to a complete snowball configuration. This result sets the ESTM limit of minimum insolation for the 
liquid-water habitability of an Earth-twin planet. 
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Fig. 9.— Variation of surface temperature as a function of latitude and orbital phase for an Earth-like 
planet with rotation rate = 0.5fl 0 (left) and il = 4^0 (right). All the remaining parameters are those 
used for the Earth reference model (Appendix B). 
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Fig. 10.— Variation of surface temperature as a function of latitude and orbital phase for an Earth-like 
planet with insolation 10% lower (left) and 10% higher (right) than the present-day value. All the remaining 
parameters are those used for the Earth reference model (Appendix B). Blank areas in the left panel indicate 
regions of the {ip, t) space with surface temperature below 200 K. 
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C.3. Atmospheric columnar mass 


In Fig. 


11 we show how T{}p,t) is affected by variations of surface pressure, the left and right panel 

Since the surface gravity is kept fixed at 


corresponding to the cases p = 0.5 bar and 4 bar, respectively. 
g = this experiment also investigates the climate impact of variations of atmospheric columnar mass, 
p/g. We find a significant difference in mean temperature and habitability between the two cases, with 
T = 274 K and h\^ = 0.62 in the low-pressure case and T = 310 K and h\^ = 1.00 in the high-pressure case. 
The mean equator-pole difference, decreases from ATep =55 K to 33 K between the two cases. The rise of 
mean temperature is due to the existence of a positive correlation between columnar mass and intensity of 
the greenhouse effect. The decrease of temperature gradient results from the correlation between p/g and 
the efficiency of the horizontal transport. Both effects have been already discussed in Paper I. Here we find 
a more moderate trend with p/g, as a result of the new formulation of D that we adopt. 

C.4. Planet radius or mass 


In Fig. 12 we show how T((p,t) is affected by variations of planet radius, the left and right panel 
corresponding to the cases R = 0.5i ?0 and 1.5i?®, respectively. In this ideal experiment we keep fixed 
the planet mean density, p = p(^, so that the planet mass and gravity scale as M (x and g (x R, 
respectively. We also keep fixed the columnar mass, p/g = {p/g)(^, by scaling p and g with R. With this 
experimental setup, the radius is the only parameter that varies in the transport coefficient D [Eqs. (231, 


(24), and (25)]. The left panel corresponds to the case M = 0.125Mq, p = 0.5bar and g = 0.5 and the 
right panel to the case M = 3.375M 0 , p = 1.5bar and g = 1 . 5 ^ 0 . 

We find that the planet cools significantly when the radius, mass and gravity increase, with a variation of 
the mean global temperature from T = 294 K to 281 K. The equator-to-pole temperature differences increases 
dramatically, from ATep =18 K to 64 K, respectively. This change of ATpp is much higher than that found 
for the same change of radius in the case of the aquaplanet (bottom-right panel in Fig. [^. The inclusion 
of the ice-albedo feedback in the present experiment amplifies variations of ATpp. In fact, the ice cover 
increases from 0% to 21% between R = 0.5 i ?0 and 1.5 R^. As a result of the variations of T and ATpp, the 
mean global habitability is significantly affected, changing from /iiw = 1-00 to 0.71 with increasing planet 
size. 


C.5. Axis obliquity 


In Fig. 


13 we show how is affected by variations of axis obliquity, the left and right panel 

and 45°, 


corresponding to the cases e = 0° and 45°, respectively. We find a modest decrease of T, from 291K to 
289K, and a significant decrease of ATpp, from 42 K to 15 K. The mean annual temperature at the poles is 
lower at e = 0° than at e = 45°, because in the first case the poles have a constant, low temperature, while 
in the latter they alternate cool and warm seasons. As a result, the ice cover decreases and the habitability 
increases in the range from e = 0° to e = 45° For the conditions considered in this experiment, the initial ice 
cover is relatively small (~ 7%) and the increase of habitability relatively modest (from 0.87 to 0.95). Larger 
variations of habitability are found starting from a higher ice cover. These results confirm the necessity of 
determining T{ip,t) and accounting for the ice-albedo feedback in order to estimate the habitability. 

At e > 45° EBM studies predict a stronger climate impact of obliquity, with possible formation of 
equatorial ice belts ( Williams fc Kasting|1997 Spiegel et al.|2009' Vladilo et al.|2013 ). The physically-based 
derivation of the coefficient D prevents using the ESTM when the equatorial-polar gradient is negative. 


because Eqs. (21) and (23) require ST > 0. In the Earth model this condition is satisfied when e < 52° 


Clearly, the climate behavior at high obliquity should be tested with 3D climate experiments (e.g. Williams 
fc Pollard|2003 Eerreira et al.|20i4 and refs, therein), being cautious with predictions obtained with EBMs. 
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Time (yr) Time (yr) 

Fig. 11.— Variation of surface temperature as a function of latitude and orbital phase for an Earth-like 
planet with surface pressure p = 0.5 bar (left) and 4bar (right). All the remaining parameters are those used 
for the Earth reference model (Appendix B). 



Time (yr) Time (yr) 

Fig. 12.— Variation of surface temperature as a function of latitude and orbital phase for an Earth-like 
planet with radius R = 0.5 i ?0 (left) and 1.5 i ?0 (right). The mean density of the planet is assumed to be 
constant, p = so that the same model predictions are also appropriate for an Earth-like planet with 

mass M = 0.125M 0 (left panel) and 3.375M 0 (right panel). The atmospheric columnar was kept constant 
by scaling p and g with R. All the remaining parameters are those used for the Earth reference model 
(Appendix B). 
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Temperature seasonal cycle (EBM oblO) 


Temperature seasonal cycle (EBM obl45) 




Time (yr) 




Time (yr) 


Fig. 13.— Variation of surface temperature as a function of latitude and orbital phase for an Earth-like 
planet with axis obliquity e = 0° (left) and 45° d (right). All the remaining parameters are those used for 
the Earth reference model (Appendix B). 



Temperature seasonal cycle (EBM PolarContInent) 




Time (yr) 


Fig. 14.— Variation of surface temperature as a function of latitude and orbital phase for an Earth-like 
planet with global ocean coverage of 0.7 and a single continent covering the remaining surface. Left panel: 
equatorial continent spread over all longitudes. Right panel: polar continent centered at latitude —90°. All 
the remaining parameters are those used for the Earth reference model (Appendix B). 
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C.6. Ocean/land distribution 

In Fig. [^we show how T((y9, t) is affected by variations of the geographical distribution of the continents. 
In these experiments we consider a single continent covering all longitudes, but located at different latitudes 
in each case. The global ocean coverage is fixed at 0.7, as in the case of the Earth. In the left panel we 
show the the case of a continent centered on the equator, while in the right panel a continent centered 
on the southern pole. The variation of is quite remarkable given the little change of mean global 

annual temperature (T = 288 K and 289 K for the equatorial and polar case, respectively). The mean annual 
habitability is almost the same in the two continental configurations (0.86 and 0.85), but in the case of the 
polar continent the fraction of habitable surface shows strong seasonal oscillations. This behaviour is due to 
the low thermal capacity of the continents and the large excursions of polar insolation. 


C.7. Long-wavelength cloud forcing 


In Fig. 15 we show how T{ip, t) is affected by variations of the long-wavelength forcing of clouds. Analysis 
of Earth data indicates a mean value 26.4W/m^ (Stephens et al. 2012), but with large excursions (e.g. 
Hartmann et aHl992 ). To illustrate the impact of this quantity on the surface temperature we have adopted 


15W/m^ (left panel) and 35W/m^ (right panel) since this range brackets most of the Earth values. The 
impact on the mean global temperature is relatively high, with a rise from T = 277 K to 295 K between the 
two cases. The ice cover correspondingly decreases from 24% to 2%, while the habitability increases from 
hiw = 0.68 to 0.95. The mean equator-pole temperature difference shows a decrease from ATep = 52K to 
36 K. This decrease of ATep is more moderate than found for variations of rotation rate, radius and axis 
tilt. 
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Fig. 15.— Variation of surface temperature as a function of latitude and orbital phase for an Earth-like 
planet with long-wavelength cloud forcing 15 W/(left) and 35 W/m^ (right). All the remaining parameters 
are those used for the Earth reference model (Appendix B). 
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Table 1 

Fiducial parameters of the ESTM 


Parameter Fiducial value 


Description 


References/Comments 


Cml50 

210 X 10® J K“^ 

Catm,o 

10.1 X 10® J 

Csolid 

1 X 10® J 

Do 

0.66 W 

Ji 

2.2 

ai 

0.18 

an 

0.70 

Ct 

-0.11 

/3 

7.98 X 10“® (°)“^ 

<OLR>ci,o 

26.4 W 

/c^ 

0.70 

fcl 

0.60 

Ac 

0.7 


Thermal inertia of the oceans^ (^2.3. 4[l 
Thermal inertia of the atmosphere‘s (q2.3.4 [ 
Thermal inertia of the solid surfacef^ 2.3.4 ' 
Coefficient of latitudinal transport ([ 2.1.1 
Modulation of latitudinal transport (^2.1.21 
Albedo oflands^ 

Albedo of frozen su rfac es and overlooking clouds 
Cloud albedo [Eq. i 28 i] 

Cloud albedo [Eq. i ^ i] 

Long wavelength forcing of clouds^ 

Cloud coverage on water 

Cloud coverage on land and frozen surface 

Ratio of moist over dry eddie transport 


Pierrehumbert 

' 2010 

Pierrehumber 

t 

2010 

Viadilo et al. 

'rr. ' 


2013[ 


Tuned to match T-latitude profile (Fi g. l^ 


Tuned to match T-latitude profile (Fig. 
Tuned to match albedo-latitude profile (.f ig. 
Tuned to match albed o-latitude pr ofile (Fig. 
Tuned'^ using Fig. 2 in|Cess i 


essl 1976 


Stephens et al. 

2 

012l 

Sanroma Pa 

le 

i2tn2 

Stubenrauch 

2013 

Sanroma 6c Pane 

2UI2 

atuDenraucn 

2UI3 

Kaspi 6c bliowman| (|2U14 

I'lg. 2) 


Representative Earth’s value that can be changed to model exoplanets with different types of surfaces or cloud properties, 
is also tuned to match the Earth’s peak of atmospheric transport at mid latitudes, 4>max (Table [^. 

‘^The parameter ct is also tuned to match the minimum value of the albedo-latitude profile. 


Table 2 

Northern hemisphere data used to calibrate the Earth model 


Quantity 

Description 

Earth value 

Model 

Units 

<T>nh 

Surface temperature 

288.61“ 

288.60 

K 

ATpE 

Pole-equator temperature difference 

40.3“ 

41.7 

K 

<^>NH 

Fraction of habitable surface 

0.851® 

0.858 


<-A>nh 

Top-of-atmosphere albedo 

0.322“ 

0.323 


<OLi?>NH 

Outgoing longwave radiation 

240.3“ 

237.6 

Wm-** 

4*111 ax 

Peak atmospheric transport at mid latitude 

5.O'* 

4.9 

PW 


^Average ERA Interim 2m temperatures | |Dee|201l[ > in the period 2001-2013. 

^Average fraction of planet surface with temperature satisfying the liquid water criterion. 
‘^Average CERES data | |Loeb et al.|2005[|200Tj > in the period 2001-2013. 

* ^Trenberth &: Caron| ( |200l| 
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